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We investigate the Taylor expansion of the baryon number susceptibility, and hence, pressure, 
in a series in the baryon chemical potential (/is) through a lattice simulation with light dynamical 
staggered quarks at a finer lattice cutoff a = 1/6T. We determine the QCD cross over coupling at 
He — 0. We find the radius of convergence of the series at various temeperatures, and bound the 
location of the QCD critical point to be T E /T c 0.94 and Hb/T < 1.8. We also investigate the 
extrapolation of various susceptibilities and linkages to finite chemical potential. 

PACS numbers: 12.38.Aw, 11.15.Ha, 05.70.Fh 



I. INTRODUCTION 



We report physics obtained in QCD with two light flavours of dynamical staggered quarks at finite temperature, 
T, and chemical potential, fi, at lattice spacing a = 1/6T. We investigate the physics at finite chemical potential 
[]| using the method of Taylor expansions that was developed in 0] and used for QCD earlier in [1, Hj] . One of the 
quantities we investigate is the radius of convergence of the Taylor series, through which we estimate the QCD critical 
point. We also investigate the dependence on /x of various other quantities of physical interest. Finally, we investigate 




FIG. 1: Series coefficients c n (L) for a quantity that diverges in the thermodynamic limit, L — > oo, at the critical point have 
the finite size behaviour shown here. The radius of convergence, R n — c„_i/c n , plateaus for n < n*(L) before rising to infinity. 
The case shown here corresponds to a singularity on the real axis, since all the c„ in the plateau are positive. 
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the linkage of quantum numbers and its dependence on T and /i. Related earlier works are [||. 

That phase transitions are rounded off by finite size effects was discovered long back by van Hove. The most familiar 
aspects are seen when simulations are directly performed in the vicinity of the critical coupling. Quantities which 
would diverge in the thermodynamic (infinite volume) limit are finite. As a result, a lattice computation never sees a 
singularity, but infers its existence from some measures. Proofs of the existence usually involve testing extrapolations: 
such are the main remaining problems at finite temperature and vanishing chemical potential. A well-developed finite 
size scaling theory can be used to study the size, L, dependence of such quantities and extract critical exponents. To 
date, the immensity of computational requirements has prevented full use of this theory for QCD. 

The study of the effect of finite size rounding of critical points on series expansions is, to the contrary, in its infancy. 
The clearest fact about such effects is the following: since quantities which should diverge in the thermodynamic limit 
merely have finite peaks at finite L, series expansions, strictly speaking, have finite radii of convergence for finite L. 
In the limit L — > oo the radius of convergence estimates the nearest critical point. The mechanism by which this limit 
is reached is straightforward. Were one to examine some estimator of the radius of convergence at order n, R n (L), 
one would find that up to some n < n*(L) the R n (L) would approach a finite value. Such behaviour was exhibited for 
the series expansion in QCD in For larger n the R n would diverge, in accordance with van Hove's theorem. With 
increasing L one would find that n„(L) approaches infinity. How the scaling of n*(L) with L codes for the critical 
exponent is at present unknown. 

Another question that can be answered by the series expansion is where the singularity lies. In the case of QCD, 
the question is whether the finite radius of convergence of the series expansion in fj, is due to a singularity at real /i. 
If it is, then the successive coefficients in the expansion are positive. At finite L one must examine the signs of the 
series coefficients for n < n*(L). If these are positive, then the singularity which limits the expansion when L — > oo 
is on the real axis. In Q this argument was used to justify the identification of the radius of convergence with the 
critical point. This argument is also implicit in Q. 

In summary, the quark number susceptibilities are Taylor coefficients in the expansion of the pressure in powers of 
the chemical potential. From the series expansion for the pressure, 

P(7> B ) = P(T) + \ X { b\t)A + ^xf{T)t& + ^x [ b\t)& + ^xi 8) (T)/i| + • ■ ■ , (I) 

we define the non-linear susceptibilities (NLS) of the n-th order, Xb ■ The second order susceptibility, also called the 
quark number susceptibility (QNS) Q, has the expansion 

X b(T,» b ) = X [ b\t) + \x ( b\t)l? b + ^x ( b\t)^ b + ^xLVVl + ■ • • • (2) 
This series is expected to diverge at the QCD critical end point. We define the radius of convergence of this series as 

,W _ 



' Xb (3) 



\n(n-l) X W 



in) 

When successive estimators for /i* are equal within statistical errors to the same value /i*, we have identified the 
plateau in the radius of convergence. This corresponds to the critical point, provided the singularity in the series 
occurs at a real value of /Lt* . In turn, this is the case when the coefficients from which the estimates are made are all 
positive. 

In the next section we present the details of the simulation and the extraction of the critical coupling. This is 
followed by a section in which we report the main results, namely, the extraction of the QNS up to the eighth order. 
This results in five terms of the series for the pressure, and four terms of the series for the baryon number susceptibility. 
Using these we report our result for the radius of convergence of the series, and extract from this our best estimate 
of the critical point. In the section after this we discuss the extrapolation of physical quantities to large chemical 
potentials. This extrapolation throws more light on the nature of the QCD critical point. 

Note that the series in eqs. (|f I2j) cannot be continued beyond /i* even when all the terms are known exactly. The 
truncated series expansion fails even faster. As a result, it becomes difficult to extrapolate physical quantities to 
large values of \j,b- One way to use the series to better advantage is well known: the method of Pade approximants. 
The existing theory of Pade approximants 0] is adapted to the case where each known series coefficient has infinite 
numerical accuracy. When coefficients are extracted through Monte Carlo estimates, and hence have statistical errors, 
new issues arise. We believe that it would be useful to extend the theory of Pade approximants in this direction. In 
the appendix we make a beginning which is adequate for the purpose of this paper. 
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TABLE I: The simulation parameters. These simulations used the R-algorithm with a step size of STmd = 0.01 and a trajectory 
length of Tmd ~ 1. For tests of accuracy, see later. 



Operator 


St = 0.01 


St = 0.001 


(P.) 
(Pt) 
(Re L) 


1.611 (1) 
1.611 (1) 
0.031 (3) 
0.293 (9) 


1.611 (2) 
1.611 (2) 
0.295 (5) 
0.291 (9) 



TABLE II: Comparison of bulk quantities, namely the spatial and temporal plaquette averages, the Wilson line and the chiral 
condensate in runs with two different MD time steps. In both cases the trajectory length was 3 MD time units and the first 
1002 MD time units were discarded for thermalization. The number of trajectories used in the comparison was 2916 for the 
larger time step and 733 for the smaller time step. 



II. SIMULATIONS 



The simulations were performed using the R-algorithm for hybrid molecular dynamics. This uses a finite step size, 
5t, for the molecular dynamics. Our main data set is generated using St = 0.01 and a total trajectory length of t = 1 
in MD time units. We performed tests of the accuracy and efficiency of the choices. 
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FIG. 2: The first figure shows \l as a function of f3. There is some finite size shift in the position of the peak at the lowest 
volumes. The second figure shows the volume dependence of the value of the peak, x™ ax as a function of N s . It is clear that 
the peak grows slower than the volume A" 3 . 
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Most of our computations were performed with St = 0.01. This was found to be adequate for computations at 
N t = 4. We checked our results at T/T c = 1.00 by running a long computation with St = 0.001 and trajectory length 
of 3 MD time units. We found complete agreement between the runs with two different time steps. In Table Hi] we 
show the comparison of bulk quantities computed in the two runs. 

Changing the trajectory length from t = 1 to t = 3 at T/T c = 0.94 ± 0.01, 1.00 ± 0.01 and 1.92 ± 0.05 did not 
change the results for thermodynamic quantities within errors. However, near T c the longer trajectories were more 
effective at reducing the autocorrelation time. For example, we found that the longest autocorrelation time at (3 C 
was Ti n t ~ 267 trajectories for T = 1, and it reduced to 36 for t = 3. As a result, the CPU time taken to produce 
decorrelated configurations is reduced by a factor of about 2.5 on taking the longer trajectories. At T/T c = 0.94±0.01 
the effective speedup, computed in the same way, was a little under a factor of two. In the high temperature phase 
the autocorrelation times were very small, and there was little to be gained by using longer molecular dynamics 
trajectories. There were no changes in thermodynamic quantities on changing the trajectory length. 

A range of (3 was scanned to locate the bare coupling at the finite temperature crossover, (5 C . The crossover was 
located by the peak of the unrenormalized Polyakov loop susceptibility, 

XL = NH(L 2 )-(L) 2 y where L = — 3 £ ReTrP(x), and P(x) = ]JU { (x,t). (4) 

Here U^(x,t) is the gauge link in the time direction at the spatial site x on the time slice t. We shall show later that 
the fourth order quark number susceptibility also peaks at the same coupling. This is closely related to the inflection 
point of the second order susceptibility which is used by various groups [9j. 

From the peak of xl we identify (3 C = 5.425 ±0.005, where the uncertainty is due to resolution, and not a statistical 
uncertainty. There is a little finite volume shift in the position of the peak of \l at the smallest volume, but no such 
shift is observed in going from N s = 18 to N s = 24 (see Figure [2]). While some volume dependence is visible in the 
peak of xl, with data from just these three volumes it is not possible to decide whether there is a crossover or a 
critical point at (3 C . However, it is possible to rule out a first order transition, since xT^f-^s definitely drops with 
increasing N s , as shown in Figure [2] 

Subsequently, T — runs were performed on lattices of size 16 4 and 24 4 on a grid of (3 to determine the scale. The 
scale determination used the value of the plaquette to obtain the renormalized gauge coupling in the MS scheme. The 
errors in the scale setting involve the uncertainty in the location of the crossover coupling, /3 C , statistical errors in 
plaquette measurements, and scheme dependence estimated by evaluating the scale also in the E and V schemes. In 
the range of temperatures within 20% of T c , the largest errors came from the uncertainty in the determination of (3 C . 
Better results can only be obtained by using larger spatial volumes. At larger temperatures, the scheme dependence 
of the scale set the largest errors. These can be reduced by going to smaller lattice spacing, i.e., to larger Nt. The 
scale setting using the crossover for Nt = 6 is compatible within errors with that obtained earlier for a similar setting 
of scales for N t — A. 



III. QUARK NUMBER SUSCEPTIBILITIES 

A quick reminder of our notation @, S| is in order. A quark number susceptibility is obtained by taking a derivative 
of the pressure with respect to the chemical potential. In two flavour QCD there are two possible chemical potentials, 
fi u and /i^- If one takes j u derivatives with respect to fi u and jd with fig, then the order of the quark number 
susceptibility is n = j u + jd- Since the u and d quarks are degenerate, and indistinguishable at fi u = (X4 = 0, we 
denote the susceptibilities by Xj u j d when j u > jd and Xjdju when jg > j u . The susceptibilities are constructed 
from expectation values of a string of 70 operators sandwiched between quark propagators. The operator O n is the 
operator with n insertions of 70 into a single fermion loop, and hence contributes only to Xno- The operators O a bc- - 
are products O a Ob and may contribute to several of the Xnm- The construction of O n on the lattice is given in 
detail in We shall discuss results for x„ m as well as the expectation values (T/V)(O a bc---) (since we discuss only 
the connected pieces of these expectation values, we have not used separate notation for that). 

The quark number susceptibilities are obtained as expectation values of fermion loops with various operator in- 
sertions @. These are evaluated as usual through stochastic estimators. The computations were optimized using 
the methods of The need to use large number of stochastic vectors has been discussed in detail elsewhere [Icj . 
We have taken 500 random vectors for each trace evaluation. With this we are able to control statistical errors on 
loops with up to six operator insertions. Even so, loops with larger numbers of insertions remain noisy. Thus, at 
T c on the largest volume, X20 gives a signal at 53<r and X40 a t 23c, whereas Xeo an d Xso gi ye signals at 5a and 3a 
respectively. For the two highest susceptibilities, this level of the signal is an improvement over the corresponding 
results with Nt = 4 at equal N s /N t . It was our experience at coarser lattice spacing that one needs lattices with 
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FIG. 3: The flavour diagonal and off-diagonal quark number susceptibilities on 6 x 24 3 (circles) and 4 x 16 4 (boxes) lattices. 

larger spatial volumes to control loops with more insertions. We see this also at the current lattice spacings, at the 
smallest volumes, even loops with six insertions are hard to control. 

In the following sections we will often compare results for N t = 4 and N t = 6. These results are meaningful only 
if they are done holding other factors fixed. We shall therefore compare the new results obtained on 6 x 24 3 lattices 
with those obtained earlier on 4 x 16 3 with the same quark mass, m/T c . 

A. Second order 

The lowest order quark number susceptibilities are shown in Figure [3J The diagonal susceptibility, X20 seems to 
show significant dependence on Nt, i.e., the lattice spacing. This is not a surprise; after all, even in the quenched 
theory a similar effect was seen 

EH- 

The off-diagonal susceptibility, xn seems to scale better with the lattice spacing. 
Note that xn takes contributions only from (T/V)(On), whereas X20 has contributions from this as well as 
(T/V)(0 2 ). The results shown in Figure [3] indicate that the quark-line disconnected operator has, at best, marginal 
lattice spacing dependence. Most of the lattice spacing dependence seen in X20 therefore comes from (T /V)(02) ■ 
This last expectation value is the response to a chemical potential on the isospin component J3 and hence was called 
X3 in some of our early papers. Both this and X20 change rapidly near T c and the "inflection point", i.e., the point 
at which the slope is maximum can be used as a corroborative measure of (3 C . Because of the numerical uncertainties 
in taking derivatives of noisy data, we will instead use the peak in the fourth order susceptibility. We discuss these 
next. 



Fourth order 



Two of the fourth-order susceptibilities are shown in Figure SJ Both X40 and X22 peak at T c . This was already 
seen in earlier simulations with Nt — 4. Within the resolution of our measurements, we see that the peak in these 
quantities comes at exactly the same coupling as the peak in \l, at both Nt — 4 and 6. Like the second order 
susceptibilities, these too have significant cutoff dependence. X3i is much smaller than either of these susceptibilities 
and shows no special structure near T c . 

The peak at T c can be resolved into a single operator expectation value, (T /V)(022) ■ This expectation value peaks 
at T c and falls off rapidly on both sides, as shown in Figure [5] Therefore it can serve as a good measure of the critical 
coupling. The expectation value (T/V){Oi), on the other hand shows a crossover near T c . One could construct yet 
another measure of the critical coupling from the point of steepest slope of this expectation value, or from its variance, 
the expectation value (T /V)(C>44) . This last quantity contributes to eighth order susceptibilities. 

Using the fourth and second order quark number susceptibilities, one can form the first two terms of the series 
expansion of X2o(A t s) El- From these coefficients can obtains the lowest order estimate of the radius of convergence 

of this series, /x* . This is shown as a function of T/T c in Figure [6l Note that the large dependence on the lattice 
spacings seen in each of the susceptibilities almost cancel out in the estimate of the radius of convergence. The radius 
of convergence has smaller dependence on the lattice spacing. 
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FIG. 4: Two of the fourth order susceptibilities on 6 x 24 3 lattices. Also shown is a comparison of X40 on 6 x 24 3 and 4 x 16 3 
lattices. 
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FIG. 5: The operator expectation values, (T/V)(e> 2 2> and (T/V){Oa) on 6 x 24 3 (circles) and 4 x 16 3 (boxes) lattices. The 
former peaks at T c whereas the latter exhibits a rapid cross over. 
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FIG. 6: The lowest order estimate for the radius of convergence of the Taylor expansion on 6 x 24 3 (circles) and 4 x 16 3 (boxes) 
lattices. 
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FIG. 7: Two of the sixth order susceptibilities on 6 x 24 3 lattices. Also shown is the expectation value which determines the 
shape of both near T c on 6 x 24 3 (circles) and 4 x 16 3 (boxes) lattices. 
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FIG. 8: The expectation value of the quark line connected operator which contributes to sixth order NLS on 6 x 24 3 (circles) 
and 4 x 16 3 (boxes) lattices. 

C. Sixth order 

The sixth order NLS is shown in Figure [7l It has been pointed out earlier that X20 has the form of a rounded 
step function, and that successively higher order NLS have the form of rounded derivatives of the step function [l2| . 
For example, the fourth order NLS has a peak. The sixth order NLS changes sign near T c and has a maximum and 
a minimum flanking the zero. This behaviour is clearly visible in Figure [7] This peculiar structure comes from the 
behaviour of the expectation value (0222), also shown in Figure [7] Note that the measurement of (O222) is noisier 
than that of (0 2 2>< 

The quark-line connected operator expectation value at this order is (Oq). This is shown in Figure [H Note that 
this has interesting structure below T c and that the structure is seen for both Nt = 4 and 6. As a result, one cannot 
use (06) or its variance for determining T c . 



D. Eighth order 



The eighth order NLS are fairly noisy below and in the vicinity of T c . This is partly due to the fact that operators 
with multiple quark loops, such as 02222, become noisier as the number of loops increases. However, single-loop 
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FIG. 9: Two of the operators which contribute to the eighth order NLS. The quark line connected operator, Og is rather noisy. 
The connected piece of the expectation value of 044 peaks at T c , as expected, and is the least noisy of all operators contributing 
at this order. 



operators such as Os are also not under sufficient control at these lattice volumes. We exhibit the expectation values 
(0$) and (©44) in Figure [H Note that there could be structure in the former below T c , but this is currently obscured 
by noise. The latter has a single sharp peak at T c , as argued before, and shows that T c identified by the peaks in xl, 
X40: X22 and (T/V)(044) are identical within the resolution of our study. 



E. Radius of convergence 



The radius of convergence of the series expansion can be used to estimate the position of the critical end point of 
QCD as before. The radius of convergence gives the distance from the origin where the expansion ceases to hold. The 
corresponding singularity lies at real ^ if the coefficients of the series expansion are all positive. As one comes down 
in T from large values of T, the series coefficients remain positive down to some lowest value of T/T c — 0.94 ± 0.01. 
At this temperature the radius of convergence is independent of the order at which it is evaluated and has value 
H* B /T = 1.8 ± 0.1. Thus, our estimate of the position of the critical end point is 

T E u E 

— = 0.94 ±0.01, and £§ = 1.8 ±0.1. (5) 

± c J- 

with a lattice spacing of 1/6T and a renormalized quark mass that corresponds to tuning m n /m p ~ 0.3, when 
the spatial size of the box is L = 4/T. In comparison, with a lattice spacing of 1/4T at the same renormalized 
quark mass and the same spatial volume, it was found that T E /T c remained unchanged whereas one had [i E ,/T E 
1 ..'i ± 0.3. Extrapolation of this result to the thermodynamic limit, L — > oo on the coarse lattice yielded an estimate 
/i^/T E = 1.1 ± 0.1, which, although statistically compatible with the finite volume result, had a lower mean. It 
would be interesting to check how much the new estimate of the QCD critical end point is lowered upon taking the 
thermodynamic limit. However, this extrapolation lies outside the scope of the current work because of the CPU 
resources needed. 



IV. PHYSICS AT FINITE fi 



A. Fluctuations and the quark number susceptibility 



Baryon number fluctuations, by an amount AB from the expectation value in a grand canonical ensemble, have a 
spectrum 

P(AB) = exp (-^^) , so ((AB) 2 )=VT XB . (6) 
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FIG. 10: Xb{hb,T e ) /(T e ) 2 obtained through various extrapolations to finite chemical potential. On the left are shown the 
series expansions to orders n — 2 and n — 4 (obtained from the 4th and 6th order NLS). The panel on the right shows the 
extrapolations using Pade approximants. 



It has been proposed that the susceptibilities be measured in event-to-event fluctuations in heavy-ion collisions 
Indeed, the divergence of the width of the spectrum of fluctuations could be one signal for the detection of the QCD 
critical point in experiments [14|. In view of this, it is interesting to estimate xb as a function of fiB along the critical 
isotherm. 

While the truncated series expansion can be used to estimate the radius of convergence, it cannot be used to 
extrapolate the susceptibility up to that point. As shown in Figure [TU1 the series expansion for xb(hb,T)/T 2 taken 
to orders n = 2 and n = 4 fail to agree long before the radius of convergence is reached; nor do they show any 
divergence at fj,*. In order to extrapolate the QNS, one must therefore find more robust techniques. 

The usual method is to convert the series to a Pade approximant (see [1 51 ] for a previous application to QCD). 
There is extensive literature on the use of these methods when the series expansion is known exactly. In Appendix 
[A] we extend this theory to the case relevant to our study, i.e., when the series coefficients are known only through 
a Monte Carlo procedure, and hence are known with certain errors. The appendix examines error propagation in 
the Pade approximants, and sets out the basic methods to control these errors. For our purposes we use the Pade 
approximants labeled Pi(i± 2 B /T 2 ) in the notation in Appendix [A] 

The Pade approximants P^(^l 2 /T 2 ) and Pl(fx 2 /T 2 ) are shown in FigurelTOl It is interesting to note that they diverge 
as Hg/T is approached. While they disagree with the series expansions as the radius of convergence is approached, 
they remain consistent with each other except very close to the divergence. Note that the errors are large near the 
divergence. This seems unavoidable, since any error in the coefficients will be magnified near the pole. There are also 
large errors beyond the pole. It should be possible to control these in future work. 

Note that in two flavour QCD one has xb = %Xbq, at all fiB, as long as the isospin chemical potential remains 
zero. So there are two independent susceptibilities, xb and XQ- I n terms of the previously computed quantities, they 
are [13], 

2 2 
Xb = g (X20 + Xn) ! and xq = (5X20 - 4xn) • (7) 

It turns out that xn remains small within errors even at larger chemical potential, so that the behaviour exhibited 
in Figure [TU] for xb is also almost quantitatively correct for xq after an overall normalization by a factor of 5/9. In 
particular, the divergence in xb is also reflected in xq- This has consequences which we deal with next. 



B. Linkage 

Earlier works have introduced quantities which measure whether two quantum numbers vary together in thermody- 
namic fluctuations [H, [13] ■ The most straightforward measure, called the linkage, utilizes diagonal and off-diagonal 
QNS in the form of the ratio 

n - r — XNM (a\ 

<~-(NM)\N — <~;(MN)\N — ) (°) 

XN 
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FIG. 11: The linkages Cnjm\v (left) and Ccbq)\q (right) as functions of T/T c at vanishing chemical potentials, as determined 
on 6 x 24 3 lattices. 
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for any two quantum numbers N and M. The linkage gives the thermal averaged amount of the quantum number 
M excited per unit N in a thermal fluctuation taking place in the grand canonical ensemble. In two-flavour QCD 
one may measure the linkage between U and D quantum numbers (conventionally +1 for quarks, -1 for antiquarks of 
the correct flavour, and zero otherwise). Also related is the linkage between the baryon number, B, and the electrical 
charge, Q. 

In Figure ITTI we show the temperature dependence of C(jjd)\u- At T = this quantity should be —2/3, since the 
lightest excitation is a pion, and the two charged pions each give a contribution of —1, whereas the neutral pion 
gives a contribution of 0. In the high-temperature phase, when the lightest excitations are quark quasi-particles, the 
linkage should vanish. We see a rapid cross over between these two regimes, with a very small but non-zero value 
being reached at T c . We also exhibit the linkage C(bq)\q- At T — this quantity is expected to vanish, since the 
lightest charged particle, the pion, has no baryon number. In an ideal quark gas, this linkage has value 1/5. One sees 
a rapid crossover between these two values in the vicinity of T c , exactly as for Cmxjjp- 

In the chiral limit, i.e., when the quark masses vanish, and a second order phase transition to occur, one would 
expect that (T/V){022) becomes infinitely peaked at T c . As a result, one expects the diagonal susceptibilities to 
become infinitely sharp, and the linkages to jump abruptly across the transition. Some part of the rounding in the 
linkages is therefore due to the fact that the quark masses are finite. However, the rounding of the crossover in 
the linkages would be a direct demonstration that there is no abrupt change from the hadronic to the quark phase: 
one may use either description over a small range of temperatures near T c . This could have implications for the 
description of hadronization in a fireball, a process which, at the moment, has a very crude description in terms of 
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the Cooper-Frye mechanism 18]. However, a part of the rounding is also due to finite volume effects, and it is hard to 
disentangle the two in our computation. It would be an interesting future computation to understand quantitatively 
what part of this slow crossover is a finite volume effect and how much is the effect of a finite quark mass. 

Since we have control over the higher order NLS, we can construct a Taylor series expansion for the linkage and 
examine its behaviour at finite chemical potential. For the analytic continuation of the linkage, we perform a jack-knife 
analysis. In each jack-knife bin the Pade approximant is evaluated at the chemical potential of interest. The mean 
of these values is used as the estimator for the continuation, and the 68% interval, evaluated non-parametrically, is 
quoted as the error bound. The results are shown in Figure fTS] for several different temperatures. 

At the highest temperature, i.e., T/T c ~ 2, the linkage Cn/D)\U is close to zero at /xb = and remains zero for 
Hb/T ~ 1. At temperatures below T c , the linkage is non-zero at (is = but evolves smoothly with the chemical 
potential. For T > T E we find a smooth increase with /xg, the linkage decreasing with larger /xg. This illustrates 
the important point that a finite radius of convergence for one susceptibility does not imply divergences in other 
quantities. 

Interestingly, the linkage C(bq)\q seems to fall marginally with increasing [is- This mild effect can be traced to the 
fact that the fourth order coefficient in the Taylor expansion of this linkage is small and negative. This results in a 
fall at large /xg. It would be useful to check whether this persists at larger volumes, and whether higher order terms 
in the expansion turn this around and cause the linkage to increase. An interesting alternative possibility is that the 
fall in C(bq)\q is physical, as is the rise in C(j/£))|£/i an d the two together imply the existence of a phase analogous to 
the quarkyonic phase at large N c [l{|. It is therefore of interest to check these results further. Unfortunately, both 
checks require massive investment of CPU resources, but are interesting enough that we hope to return to this soon. 



V. SUMMARY 



We have examined QCD with two flavours of dynamical staggered quarks at finite temperature with lattice spacing 
a = 1/6T and bare quark masses tuned to m/T c = 0.1. This quark mass is expected to correspond to m^jmp = 0.3, 
and hence our new results are directly comparable to the older results which were obtained on a coarser lattice with 
a = 1/4T [3j. Our simulations were performed on lattices with size LT < 4, where L is the spatial size of the box. 
We used the R-algorithm with a step size in MD time units of 0.01. We have checked that decreasing this by an order 
of magnitude to 0.001 does not change thermodynamic results (see Table ITT|) . Similarly, we have checked that the 
physics results remain unchanged when the trajectory length is changed. 

We identified the cross over at vanishing chemical potential through the Polyakov loop susceptibility, \l, (see Figure 
[5]) and then cross checked this through two measures related to the QNS. One is the peaking of X40 and X22 (see Figure 
HJ) which is related to the "inflection point" of the QNS. The other is the peaking of the operator (T/V){04a), (see 
Figure [9]) which is related to a similar inflection point in (T/V){0^). These measures are consistent with each other 
within the accuracy of our computations. The scale setting using this identification of the cross over is consistent 
with the earlier scale setting using coarser lattice spacing Q. 

We presented results for the NLS up to eighth order. There are clear lattice spacing effects, as expected. These 
are roughly consistent with earlier determinations of some of these quantities in quenched theory. While the lattice 
spacing artifacts for the NLS are very large, sometimes as much as 100%, the effect on the radius of convergence is 
much smaller (see Figure [5]). Our estimate of the critical point of QCD is based on this radius of convergence. The 
critical point occurs when the radius of convergence identifies a singularity on the real axis, through the fact that the 
series coefficients are all positive. Caveats on this are presented in the introduction. Our estimate of the critical point 
using finite volume data is (see eq. [5J 

T E u E 

= 0.94 ±0.01, and ^ = 1.8 ±0.1. 

T c ' T E 

This should be compared with our earlier estimate on the same lattice volume and same renormalized quark mass 
which gave fi§/T E = 1.3 ±0.3. This is a change of about 26%, and is statistically significant. Extension of our results 
to larger volumes is outside the scope of this work. In simulations with a = 1/4T, the estimate of /x|j dropped by 
about 16% on extrapolating to infinite volume. 

The series expansion is a good tool for extracting the radius of convergence, and, through it, the critical point. 
However, as we show in Figure ITOl it is a bad tool to extrapolate physical results to high zx^/T. Pade approximants 
adjusted to give the same series expansion seem to perform better, even after taking into account the propagation 
of statistical errors. One sees the divergence of the susceptibility at the critical end point, something that the series 
expansion misses altogether. 

We also examined the linkages C(jjb,)\u and C(bq)\q (see Figure [TT1) . At vanishing chemical potential they show 
a rapid cross over from the values expected in the hadronic phase to those expected for the nearly ideal quarks. 
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The rounding of this transition is closely related to the question of how sharp the hadronization transition can be in 
heavy- ion collisions. A discussion of the issue was presented in Section HVB I 

The measurements of the linkages were extended to finite chemical potential using Pade approximants. Unlike the 
QNS, they evolve smoothly through the critical point. Interestingly, on isotherms below T Cl with increasing the 
linkage between U and D quantum numbers changes towards the ideal quark gas, whereas the linkage between B and 
Q changes away from the quark gas. This could indicate the presence of a quarkyonic phase of QCD matter, although 
technicalities need to be sorted out before one can establish this. 

The computations were performed over the last two years on the Cray XI of the Indian Lattice Gauge Theory 
Initiative (ILGTI) at TIFR. We thank Ajay Salve for single-handedly taking care of the machine during this extended 
period. 



APPENDIX A: PADE APPROXIMANTS 



We follow Baker's definition [8| of a Pade approximant [2fJ,|21|. The series expansion 

In{x) =c +c lX + --- c n x n + 0(x N+1 ) 
evaluated to order N can be used to define the Pade approximant of order L/M, 



PmUn(x)] = -f 



B M (0) = 1, Btr(x)f N (x) - AUx) = 0{x L+M+ ") 



(Al) 



(A2) 



where Aj^(x) and B M (x) are polynomials in x or order up to L and M respectively. From the matching condition it 
follows that L + M < N. Introduce the notation — 



A M (x) = a + a x x + • ■ ■ + a L x L : B M {x) = 1 + b%x + • ■ ■ + b M x 



A I 



(A3) 



Then, writing out the matching condition order by order, one obtains the Pade approximants by solving first for the 
denominator 



f CL-M+l CL-M+2 ■ ■ ■ CL \ / &M \ 



CL-M+2 CL-M+3 



cl+i 



\ 



JM-1 



(-L 



CL+1 



cl+m- 



J \ h J 



I CL+1 \ 
CL+2 



V C L +M I 



(A4) 



(with the convention that Cj = for j < 0) and then for the numerator 



a = c , 

ai = ci + &ic , 



(IL 



min(L,M) 



biC L - 



(A5) 



The practical importance of Pade approximants arises from the fact that if the series /at has a radius of convergence 
R as N — *■ oo, then the series expansion is reliable only for x < R, whereas the Pade approximants can be used for 
analytic continuation beyond this. Much of the standard theory of Pade approximants deals with the cases when the 
coefficient matrix in eq. (|A4j) has vanishing determinant, and the information which can then be extracted. 

Here we concentrate on a different problem — that of controlling errors when the series coefficients are obtained 
by a Monte Carlo program, and hence have a given statistical distribution. We found no discussion of this in the 
literature, although it is likely that sporadic attempts to answer related questions have been made in the past. These 
questions become important now that new developments in QCD at finite chemical potential lead us to analyze series 
coefficients obtained in a Monte Carlo process. 

When the Pade coefficients are well-defined, the joint probability distribution of the series coefficients can be 
transformed into that of the coefficients of the Pade approximant using the usual Jacobian formula — 



T ? M( a o, ai, • • • , a,L,bi,b 2 , b M ) = V(c , c\, ■ ■ ■ , c l +m)J, 



where 



d(a ,ai, 



,a L ,bi,b 2 ,- ■ ■ ,b 



M) 



9{cq,ci, • • • , cl+m) 



(A6) 
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Take the example of Pj 3 , where eto = Co and b\ = Ci/co, so that J = ao- Assume that cq and ci are drawn from 
independent Gaussian distributions of unit mean — 



Pi°(c ,ci) 



1 



27T(7o cr l 



exp 



i r (c 

2 



I) 2 , 



2a? 



(A7) 



Then the joint distribution of the Pade coefficients can be written down. The marginal distribution of the Pade 
coefficient bi, being the ratio of two Gaussian distributed numbers, is well known [13] and given by 



where 



27rA(fei) 



27T(ToC r l^ 2 (^l) 



(A8) 



A 



















4) 


M = 













A 2 
4 ' 



Note that A(±oo) = 1/crf, and hence /i(±oo) is a finite number which depends only on <xo,i- As a result, the marginal 
distribution of b\ is not exponentially damped at infinity. The power-law damping comes from the factor 1/v 2 ~ 1/fr 2 - 
Clearly this distribution has a well-defined expectation value for b\, but the variance and higher cumulants do not 
exist. Thus, statistical measurements of b\ are not subject to the central limit theorem. 

A similar phenomenon occurs with any Py . Assume that the series coefficients are statistically independent and 
drawn from a Gaussian of unit mean, then the probability distribution of the Pade coefficients can be written as 

/L + l 



Vi(a ,ai,---,a L ,bi) 



n 

,i=0 



1 



1 I Q 

J exp — — 



(A9) 



where J is the Jacobian of the transformation given by Cj = Ylj<i a i-j(~ ^lY for i < L and c^+i = c^bi, and Q is 
quadratic form obtained by transforming the arguments of the Gaussians. 
The Jacobian of this transformation is 



Jl = 



1 

-bi 

(-M 2 

(-h) L 
bi(-h) L 




1 



{-bi) L - 1 
bii-h)^ 1 





1 

(-bi) L - 2 
bi(-h) L - 2 



c' 2 {bi) 



c' L (h) 



(A10) 



where c[{bi) is the derivative of Ci with respect to b\. Multiplying the second row from the bottom by b\ and 
subtracting that from the last row gives 



Jl 



1 

-h 
(~bi) 2 

(-bi) L 






1 

-bi 







(-bi) L - 2 






c[ 
c' 



c L 



CL 



(All) 



i=l 



where we have used the relation cl+i = b\c^ to write c' L+1 — b\c' L = cl- 

The quadratic form in the argument of the exponent can be manipulated into a particularly useful form by com- 
pleting the squares — 



Q ee J2 ~ 2 =a T Qa + 2b T a+^4 = ( a -") T Q( a - s ) 



L+l 



/i. 



L+l 



where 



/' : 



5^ rr? 



a Q& and a = Q 1 b, 



(A12) 



i=l 



where the real symmetric matrix Q and the vector b can be easily written down. We do this next for the special case 
when all the Oi are equal to a. 
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Define the sequence of polynomials 

Pj (bx) = 1 + 6? + 6* + ••• + 6^' = l + b 2 lPj ^(h), 



In terms of these, one writes 



Q 



qj (h) = 1 - &! + b\ - ■■■ + (-fix)* = 1 - hq^h). 



/ Pl+i -&iPx 6?Pl-i ■ ■ \ 



(A13) 
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(A14) 



In order to find the determinant of Q, we do the following row operations — starting from the top, add to each row 
b\ times the next row. This reduces the determinant to a lower triangular form 



Det Q = \ 
u 



1 

-h 
b\ 



o 
l 



Pi{h 



(A15) 



(-60 L pi (-60 i - 1 Pi {-b!) L - 2 Pi ■■■ pi 
The solution of the equation Qa = b can be obtained by the same operations. They yield the reduced equation 
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(A16) 



where we have used the relation qi(x) = 1 — xqi-\(x), to reduce the vector b. This gives- 

1 + h 



1 + bi 



(A17) 



Finally, 



L + l , T _ L + l 
— ^ b a 



4- - y>*(&o + 6i?i-i(60] - T-rn- 



(A18) 



Clearly, remains finite in the limit b± — » ±oo, so that exp(— /z/2) does not damp the marginal distribution of b±. 
In fact, that damping comes from the factor of 1/DetQ = 1/(1 + 6f). As a result, b\ has well defined mean but its 
variance is undetermined. Thus, estimators of b\ evade the central limit theorem. 

Nevertheless, the situation is pretty well under control. The appropriate question to ask of a distribution such as 
that in eq. (|A9[) in the context of parameter estimation is not the value of the variance, but an appropriate measure 
of the variation in the estimate. One could quote the width at half maximum, or the limits such that 68% of the 
probability lies within these limits. Along with this one asks, if we make ./V measurements of b\ then how does such 
a measure of variation change with N. 

A numerical investigation shows that when <7o = a% = 1, the modal value is b = 0.345897. Since the distribution is 
skew, the modal value and the mean are different. The full width at half maximum is contained in —0.3485 < b\ < 
1.26641, and this range contains 56.1% of the integral. The 68% probability interval is —0.575 < bi < 1.38. Either of 
these ranges can be quoted as an estimate of the error in the modal value. 

To answer the question about the distribution of means, we use the characteristic function. If f(x) is the distribution 
of x, then the Fourier transform f(x) is called the characteristic function. Since f(x) is non-negative and integrable, 
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being a probability distribution, it is also square integrable, so that the characteristic function exists. The characteristic 
function of the mean of N numbers, fipf, is 



Xn{w) = \ dfiN exp(zw^jv) 



N 



dx i f{x l ) 



N 



S I dxi — Nun 



f N (- 



(A19) 



Fourier transforming this gives the distribution of fi^. While this general method remains valid for the distributions 
Vi , above, it does not seem possible to perform the Fourier transformations in closed form. So, instead of writing 
down an impenetrable formula for the distribution of means, we investigate useful subsidiary questions. 
Define the skew of a distribution by 



1 



(A20) 



where x m denotes the modal (most probable) value, and (x) is the mean. The skew is nonzero for every skewed 
distribution, being positive if the distribution is skewed to the left and negative otherwise. For the distribution V®, 
we found S = 1.3 • ■ ■. For the distribution of means of N values, S decreases. A Monte Carlo estimate for ao = cr± = 1 
indicates that S « 3.5/yN. This estimate was obtained using values of TV between 1 and 100. A similar result was 
obtained for a measure of skewness that compares the median and the mean in a manner analogous to eq. (|A20[) . 

At the median of a distribution, the cumulative distribution becomes equal to 0.5. The errors on the median can 
be defined as the points at which the cumulative distribution is either 0.34 above or below. The distribution of the 
means of N numbers narrows rapidly, and we find in a Monte Carlo estimate that both these intervals decrease as 

i/Vn. 

In conclusion, for the estimation of b\ and confidence limits on the estimate, it matters little that the central limit 
theorem does not hold. The mean is well defined, and its difference with the mode and median scale with a factor 
of 1/yN. The 68% confidence limits also scale as l/\fN. We therefore quote the mean and 68% confidence limits 
on it as estimators for the Pade coefficients. These estimators are easy to incorporate into jack-knife and bootstrap 
analyses. In parts of our analysis the estimators of the series coefficients are also not Gaussian distributed; even so, 
the non-parametric statistical analysis outlined here suffices. 
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